﻿using System;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using NUnit.Framework;

namespace dnAnalytics.Tests.LinearAlgebra.Decomposition
{
    public abstract class HouseholderTests
    {
        private const string matrixFile = "./TestData/Matrices/random_real_general_array_100.mtx";
        private const string singularFile = "./TestData/Matrices/gear_integer_general_coordinate_100.mtx";
        private const string tallFile = "./TestData/Matrices/random_real_general_array_20_10.mtx";
        protected abstract Matrix GetMatrix(string file);
        protected abstract Matrix GetMatrix(int rows, int cols);
        protected abstract Vector GetVector(int count);

        [Test]
        [ExpectedException(typeof (ArgumentNullException))]
        public void ConstructorNull()
        {
            new Householder(null);
        }

        [Test]
        [ExpectedException(typeof (InvalidMatrixOperationException))]
        public void WideMatrixThrowsInvalidMatrixOperationException()
        {
            new Householder(new DenseMatrix(3, 4));
        }

        [Test]
        public void NotFullRank()
        {
            Matrix singularMatrix = GetMatrix(singularFile);
            Householder qr = new Householder(singularMatrix);
            Assert.IsFalse(qr.IsFullRank());
        }

        [Test]
        public void FullRank()
        {
            Matrix matrix = GetMatrix(matrixFile);
            Householder qr = new Householder(matrix);
            Assert.IsTrue(qr.IsFullRank());
        }

        [Test]
        public void SquareFactors()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);

            Matrix q = qr.Q();
            Matrix r = qr.R();

            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
            q.Multiply(r, result);
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    double error;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - squareMatrix[i, j])/squareMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }
        }

        [Test]
        public void TallFactors()
        {
            Matrix tallMatrix = GetMatrix(tallFile);
            Householder qr = new Householder(tallMatrix);

            Matrix q = qr.Q();
            Matrix r = qr.R();

            Matrix result = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
            q.Multiply(r, result);

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    double error;
                    if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - tallMatrix[i, j])/tallMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.AcceptableError);
                }
            }
        }

        [Test]
        public void Determinant()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);

            double expected = 5.006656015927541e248;
            double actual = qr.Determinant();
            double error = System.Math.Abs((actual - expected)/expected);
            Assert.IsTrue(error < Constants.AcceptableError);
        }

        [Test]
        public void SingularDeterminant()
        {
            Matrix singularMatrix = GetMatrix(singularFile);
            Householder qr = new Householder(singularMatrix);
            Assert.AreEqual(0, qr.Determinant(), Constants.AcceptableError);
        }

        [Test]
        public void Solve_SquareMatrix()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Matrix identity = GetMatrix(squareMatrix.Rows, squareMatrix.Rows);
            for (int i = 0; i < identity.Rows; i++)
            {
                identity[i, i] = 1.0;
            }

            Householder qr = new Householder(squareMatrix);
            Matrix result = qr.Solve(identity);

            result = squareMatrix*result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.AcceptableError);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        public void Solve_ArgumentNullExcpetion()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);
            Assert.Throws<ArgumentNullException>( ()=> qr.Solve((Matrix)null));
            Assert.Throws<ArgumentNullException>(() => qr.Solve(squareMatrix, null));
            Assert.Throws<ArgumentNullException>(() => qr.Solve(null,squareMatrix));

        }
        [Test]
        public void Solve_NotConformable()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);
            Matrix matrix = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
            Assert.Throws<NotConformableException>(() => qr.Solve(matrix));
            Assert.Throws<NotConformableException>(() => qr.Solve(squareMatrix, matrix));
            Assert.Throws<NotConformableException>(() => qr.Solve(matrix, squareMatrix));
            Matrix other = GetMatrix(squareMatrix.Rows, squareMatrix.Rows+1);
            Assert.Throws<NotConformableException>(() => qr.Solve(other, squareMatrix));
        }

        [Test]
        public void SolveVector()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);   
            Householder qr = new Householder(squareMatrix);
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0)/(i + 2.0);
            }

            Vector result = qr.Solve(vector);
            result = squareMatrix*result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0)/(i + 2.0), result[i], Constants.AcceptableError);
            }
        }
        
        [Test]
        public void SolveVector_ArgumentNullExcpetion()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);
            Vector vector = GetVector(squareMatrix.Rows);
            Assert.Throws<ArgumentNullException>(() => qr.Solve((Vector)null));
            Assert.Throws<ArgumentNullException>(() => qr.Solve(vector, null));
            Assert.Throws<ArgumentNullException>(() => qr.Solve(null, vector));

        }
        
        [Test]
        public void SolveVector_NotConformable()
        {
            Matrix squareMatrix = GetMatrix(matrixFile);
            Householder qr = new Householder(squareMatrix);
            Vector vector = GetVector(squareMatrix.Rows);
            Vector vectorP1 = GetVector(squareMatrix.Rows + 1);
            Assert.Throws<NotConformableException>(() => qr.Solve(vectorP1));
            Assert.Throws<NotConformableException>(() => qr.Solve(vector, vectorP1));
            Assert.Throws<NotConformableException>(() => qr.Solve(vectorP1, vector));
        }
    }
}